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Perturbation theory for bright spinor Bose-Einstein condensate sohtons 
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We develop a perturbation theory for bright solitons of the F — 1 integrable spinor Bose-Einstein 
condensate (BEC) model. The formalism is based on using the Riemann-Hilbert problem and 
provides the means to analytically calculate evolution of the soliton parameters. Both rank-one 
and rank-two soliton solutions of the model are obtained. We prove equivalence of the rank-one 
soliton and the ferromagnetic rank-two soliton. Taking into account a splitting of a perturbed 
polar rank-two soliton into two ferromagnetic solitons, it is sufficient to elaborate a perturbation 
theory for the rank-one solitons only. Treating a small deviation from the integrability condition 
as a perturbation, we describe the spinor BEC soliton dynamics in the adiabatic approximation. 
It is shown that the soliton is quite robust against such a perturbation and preserves its velocity, 
amplitude, and population of different spin components, only the soliton frequency acquires a small 
shift. Results of numerical simulations agree well with the analytical predictions, demonstrating 
only slight soliton profile deformation. 

PACS numbers: 03.75.Lm, 03.75.Mn, 02.30.Ik 
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I. INTRODUCTION 

Bright and dark solitons in quasi-one-diniensional 
Bose-Einstein condensates (BECs) , observed experimen- 
tally [li 0, S S] ' ^^"6 expected to be important for various 
applications in atom optics [5|, including atom interfer- 
ometry, atom lasers, and coherent atom transport. Re- 
cent experimental and theoretical advances in BEC soli- 
ton dynamics are reviewed in Refs. [g,0, Q- 

Spinor BEC of alkali atoms [3, \l^ with a purely opti- 
cal confinement, along with the two-component conden- 
sate [nl, [13, [13 J represents an example of the condensate 
with internal degrees of freedom which endow the soli- 
tons with vectorial properties. Modulational instab ility 
in the spinor BEC model was investigated in Ref. [14| . 
and some exact solutions and their stability were stud- 
ied in Ref. [15|. Vector gap solitons and self-trapped 
waves were identified in the spinor BEC model loaded 
into one-dimensional optical lattice potential [la|. Re- 
cently bright-dark soliton complexes in this model have 
been found [17| by reducing it to the completely inte- 
grable Yajima-Oikawa system |18| . 

Wadati and co-workers found [l^ that the three- 
component nonlinear equations describing the BEC with 
the hyperfine spin F = 1 admit the reduction to another 
integrable model - the 2x2 matrix nonlinear Schrodinger 
(NLS) equation, after imposing a constraint on the con- 
densate parameters. Both bright ^22,] and dark [23| soli- 
tons possessing properties of true solitons of integrable 
equations have been found. The formalism of the inverse 
scattering transform for the matrix NLS equation under 
non- vanishing boundary conditions was developed in Ref. 



•Electronic address: 'doktorov@dragon.ba s-net.by| 
^Electronic address: jwang@cems.uvm.edul 
■••Electronic address: |jyang@cems. uvm.edu 



r 



[2J| and extended in Ref. ^23| to describe bright spinor 
BEC soliton dynamics on a finite background. The full- 
time description of the modulational instability develop- 
ment in the integrable spinor BEC model was given both 
numerically [15] and analytically [2f|. 

Integrable models provide a very useful proving ground 
for testing new analytical and numerical approaches to 
study such a complicated system as the spinor BEC. At 
the same time, the integrability conditions impose spe- 
cific restrictions on the parameters of the model which 
can conflict with actual experimental settings, despite the 
fact that the effective interaction between atoms in BEC 
can be tuned, to some extent, by the optically induced 
Feshbach resonance [20, l21[. Besides, in experiment it is 
impossible to exactly hold the conditions between param- 
eters which assure integrability of the model. Therefore, 
sufficiently general analytical results concerning the full 
(nonintegrable) model with realistic parameters would be 
of importance. 

As a step in this direction, in the present paper we de- 
velop a perturbation theory for the integrable spinor BEC 
model. Evidently, small disturbance of the integrability 
condition can be considered as a perturbation of the inte- 
grable model. Our formalism is based on the Riemann- 
Hilbert (RH) problem associated with the spinor BEC 
model. The main advantage of the proposed method is 
its algebraic nature, as distinct from the method using 
the Gel'fand-Levitan integral equations [23]. The ap- 
plication of the RH problem for treatin g pe rturbed soli- 
ton dynamics goes back to Refs. [2^, ^^. The mod- 
ern version of the perturbation theory in terms of the 
RH problem has been developed in a series of papers 
[30 , [3 l| . [3 3. \33i . [34| , with its most general formulation in 
Ref. [3a ]. Another version of the soliton perturbation 
theory (the direct perturbation theory) has been devel- 
oped on the basis of expanding perturbed solutions into 
sqtiared eigenfunctions of the linearized soliton equations 
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As was shown by Wadati and co-workers [T^, l23 |. 
bright soHtons in the integrable spinor BEC model can 
exist in two spin states - ferromagnetic (non-zero total 
spin) and polar (zero total spin). Energy of the polar 
soliton is greater than that of the ferromagnetic soliton. 
Moreover, the polar soliton demonstrates a two- humped 
profile in a wide range of its parameters. Our numerical 
simulations revealed that the polar soliton is unstable un- 
der the action of a perturbation of a rather general form 
and splits into two ferromagnetic solitons. This fact is 
crucial for the development of a perturbation theory for 
spinor BEC solitons. 

The paper is organized as follows. After formulating 
the model in Sec. |TT1 we introduce in Sec. IIIII analytic 
solutions of the associated spectral problem, in order 
to formulate in Sec. IIVI the RH problem. Solving this 
problem, we derive in Sec. |V] bright soliton solutions of 
the integrable spinor BEC model, both for the rank-one 
and rank-two projectors. The rank-one soliton is char- 
acterized by the familiar hyperbolic secant profile, while 
the rank- two soliton has a more complicated form [l9|. 
Two types of the rank-two solutions are exactly ferro- 
magnetic and polar solitons. We prove that the ferro- 
magnetic rank-two soliton is equivalent to the rank-one 
soliton. In virtue of the fact that the perturbed polar 
soliton splits into two rank- two ferromagnetic solitons, 
it is sufficient to develop a perturbation theory for the 
rank-one soliton only. This is performed in Sec. IVII 
We derive evolution equations for the soliton parameters 
which exactly account for the perturbation and serve as 
the generating equations for iterations. Section |VlT] con- 
tains a description of the soliton dynamics in the adia- 
batic approximation of the perturbation theory. We show 
analytically that a ferromagnetic soliton is quite robust 
against a small disturbance of the integrability condi- 
tion, the only manifestation of the perturbation action is 
a minor shift of the soliton frequency. Numerical simula- 
tions of the perturbed spinor BEC equations are in close 
agreement with the analytical predictions revealing only 
a small soliton shape distortion and little perturbation- 
induced radiation. Section FVlIII concludes the paper. 



coupled Gross-Pitaevskii equations [2^, l3£ 
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(2.1) 



co|$oP*o + 2c2$+$-$5, 



where the constant parameters cq — {go + 2g2)/3 and 
C2 = (32 — 5o)/3 control the spin-independent and spin- 
dependent interaction, respectively. The coupling con- 
stant gf (/ = 0, 2) is given in terms of the s-wave scat- 
tering length af in the channel with the total hyperfine 
spin /, 



9f 



a/ 



ma : 



a±_ 



Here a± is the size of the transverse ground state, m is 
the atom mass, and C = —^(1/2) « 1.46. 

It was noted in [1^ that Eqs. (|2.ip are reduced to an 
integrable system under the constraint 



Co = C2 



-c< 0. 



(2.2) 



The negative C2 means that we consider the ferromag- 
netic ground state of the spinor BEC with attractive in- 
teractions. The condition (|2.2p . being written in terms 
of gf as 2go = —.92 > 0, imposes a constraint on the 
scattering lengths: a± — 3Caoa2/(2ao -I- 02). Redefin- 
ing the function $ as $ — > {(j)+, V24>o, 4'-)'^ : normalizing 
the coordinates as i -^ {c/h)t and x -^ {y^2mc/h)x, and 
accounting for the constraint (|2.2p . we obtain a reduced 
system of equations in a dimensionless form: 

idt(j)±+dlcj)±+2 (|0±|2 + 2|0on </'±+2(/);0g - , (2.3) 



tdt(t>o+dl4o+2 (|0+|2 + \M^ + l^-P) <Po+2<P+<P*o<t>- = . 

After arranging the components (j)± and 0o into a 2 x 2 
matrix Q, 



II. MODEL 



We consider an effective one-dimensional BEC trapped 
in a pencil-shaped region elongated in the x direction 
and tightly confined in the transversal directions. The 
assembly of atoms in the hyperfine spin F = I state 
is described by a vector order parameter ^{x,t) — 
{^+{x,t),^oix,t),^-{x,t))'^ , where its components cor- 
respond to three values of the spin projection m,F = 
1,0,-1. The functions ^± and $0 obey a system of 



Q H :; ;" ) , (2.4) 



we transform Eqs. (|2.3p to the integrable matrix NLS 
equation 



idtQ + dlQ + 2gQtg ^ 



(2.5) 



The matrix NLS equation (|2.5p appears as a compatibil- 
ity condition of the system of linear equations j27| 






(2.6) 
(2.7) 



where A = diag(— 1, —1, 1, 1) 
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,V^2kQ 
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(2.8) 

and fc is a spectral parameter. Equation ()2.6p (the spec- 
tral problem) enables us to determine initial spectral data 
from the known potential Qq-, while Eq. (j2.7p governs the 
temporal evolution of the spectral data. A new solution 
of Eq. ([23]) [and hence of the BEC equations ([2?3|) ] is ob- 
tained as a result of the reconstruction of the potential 
Q from the time-dependent spectral data. 



III. JOST AND ANALYTIC SOLUTIONS 

To determine the spectral data, we introduce matrix 
Jost solutions J±{x,k) of the spectral problem (|2.6p by 
means of the asymptotes Jj- — > 1 as a; — *■ ±cxd. Since 
trA = 0, we have det J± ~ 1 for all t. Being solutions 
of the first-order equation (|2.6p . the Jost functions are 
not independent but are interconnected by the scattering 
matrix S: 

J_ = J+ESE-\ E = cxp{ikAx), det S = 1. 

Besides, the Jost solutions and the scattering matrix obey 
the involution property. Indeed, since the potential Q is 
anti-Hermitian, we obtain 



Jl(r) = j±i(fc). 

Similarly for the scattering matrix: 



(3.2) 
(3.3) 



Note that the scattering matrix is defined for real k. 

For the subsequent analysis, analytic properties of the 
Jost solutions are of primary importance. Let us repre- 
sent the matrix Jost solution J as a collection of columns: 
J = (jl^J, jl^l, j[^l, jW), and consider the first column. 
Rewriting the spectral equation (|2.6p with the corre- 
sponding boundary conditions in the form of the Volterra 
integral equations, we obtain a closed system of equations 
for entries of the first column: 



1+ / dx'{(t)+J^3i + (j)oJ-ii){x'), 
da;'((/)o J-31 + (f)-J^4i){x'), 
" dx'(0;j_n + 0S^-2i)(x')e'*'=("-"'), 

OO 

' dx'(0S J-ii + ct^*_J-2i){x')e^'''^^-^'\ 



J-ii 

J-21 
J-31 

J-41 



The last two integrands point out that the column J_ is 
analytic in the upper half-plane C+, where Imfc > 0, and 
continuous on the real axis Imk = 0. This can be proved 



in the same way as for the scalar NLS equation, under 
the condition of sufficiently fast decrease of the potential 

[21 

Q at infinity. Similarly we obtain that the column Ji 
is analytic in C-|- as well, while the two other columns 
jI and J_ are analytic in the lower half-plane C_ and 
continuous on the real axis Imk — 0. As regards the 
matrix solution J_|_, its first and second columns J^ and 

[21 

— — ^--^'' '-- '^ while the third and forth ones 

. Therefore, the matrix 



J_^ are analytic in 

J^^l and 
function 



j\ and j\ are analytic in 



V'+ = (j™,ji'Ui'Ui'') (3.4) 



solves the spectral equation (|2.6p and is analytic as a 
whole in C+. 

It is not difficult to see from Eqs. (|3.ip and (|3.4p that 
the analytic solution -0+ can be expressed in terms of the 
Jost functions and some entries of the scattering matrix: 



^-4 



J+ES+E-^ = J^ES-E- 



(3.5) 



where 



S+{k) 



SU Si2 0' 

S21 S22 

S31 S32 1 

S41 S42 1 



S-{k) 



1 si, si. 



'32 ^42 
33 *43 



1 
S 



U L) s^^ S44 



(3-6) 
In writing the expression for 5*- we use the involution 
(|3.3p . These upper and lower block-triangular matrices 
S± factorize the scattering matrix lS9\: SS^ — S+. Be- 
sides, it follows from Eq. p.Sp and det J± — 1 that 



(3.7) 



detip+ = iTi\ = m_ 



where rrv^ (m_ ) is the second-order principal upper 
(lower) minor of the scattering matrix. 

To obtain the analytic counterpart of ip^ in C_, we 
consider the adjoint spectral equation 



d^K±^ik[k,K±]-K±Q 



(3.8) 



with the asymptotic conditions K± — > 1 at a; ^ ±00. 
The inverse matrix J^^ can serve as a solution of the 
adjoint equation (|3.8p . Now we write a closed system 
of integral equations for rows of the matrices K±. For 
example, the first row ^-[11 obeys the equations 



K^ 



11 



= 1+/ dx'(0;X_i3+0S^_14)(x'), 



X_i2 = / dx'(0SX_i3 + 0*_X-i4)(a:'), 

J —00 

i^-13 = - r da;'(0+X_n + (/)oif-i2)(a;')e-'^'^'"-"'\ 

J —00 

i^-14 = - r da;'(<^oi^_n + 0_X_i2)(x')e-2»'=(^-^'). 



It is seen that the row i^-fii is analytic in ' 



Similarly, 



the second row ^-[2] is analytic in C_, too, and the 
rows if- [3] and -ft^-[4] are analytic in C+. For the matrix 
solution A% we find that the rows ^+[1] and ^+[2] ai'e 
analytic in C+, while A'^ra and A'^mi are analytic in C_. 
Therefore, the matrix function 



V--^ = {K-[ih K_[2], K+[3], K+[i{)' 



(3.9) 



solves the adjoint equation p.Sp and is analytic as a 
whole in C_. Similar to ip+, the function i/jZ is ex- 
pressed in terms of the Jost solutions and the scattering 
matrix: 



^-1 = ET+E-^J^^ = ET.E-^JZ^, 
where the matrices T±, 



(3.10) 



^11 ■'21 ^31 ''41 
rp _ I •S12 ■S22 *32 '^42 
+ '0010 
1 



T_ = 



10 

10 

S31 S32 S33 S34 

S41 S42 S43 S44 



provide one more factorization of the scattering matrix: 
T- = T+S. As in Eq. dSj]), we can write 



detV^. 



m 



(2), 



,(2) 



(3.11) 



Note that the analytic solutions satisfy the involution 
property as well: 



^;(fc) = vi'(fc*)- 



(3.12) 



This property can be taken as a definition of the analytic 
function ?/;Z from the known analytic function ip^. 



IV. THE RIEMANN-HILBERT PROBLEM 

Hence, we constructed two matrix functions ■(/'+ and 
ipZ which are analytic in complementary domains of the 
complex plane and conjugate on the real line. Indeed, it 
follows from Eqs. p.Sp and (|3.10p that ijj± obey the 
relation 



iljZ\k)iP+{k) =EG{k)E- 



Imfc = 0, 



where 



Gr — 1+0+ — l—O — 



1 





«31 


si. 
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*32 


S42 


S31 


S32 


1 





S41 


S42 





1 



(4.1) 



(4.2) 



Equation (|4.ip determines a matrix Riemann-Hilbert 
problem, i.e. a problem of the analytic factorization of a 
nondegenerate matrix G in (j4.2p . given on the real line, 
into a product of two matrices which are analytic in com- 
plementary domains C±. The RH problem (|4.ip needs a 
normalization condition, which is usually taken as 



ip±{x, fc) — > 1 at |fc| — > 00. 



(4.3) 



The analytic matrix functions ip± can be treated as 
a result of a nonlinear mapping between the potential 
Q{x) and a set of the spectral data which uniquely char- 
acterizes a solution of the RH problem (|4.ip and (|4.3p . 
Conversely, the potential can be reconstructed from an 
asymptotic expansion of i^±{x,k) for large k. Indeed, 
writing ip± as 



ip+{x,k) = t + k 



-VA^^ 



0{k- 



■ipZ\x, k) = t + fc-Vi^^ + 0{k-^) 

and inserting these expansions into Eqs. (|2.6p and 
we obtain 



Q 



/,(i) 



i[A,iP'+>]=i[A,iP'_ 



(i)i 



(4.4) 



Hence, having solved the RH problem, we can find solu- 
tions of the BEC equations. 

In general, the matrices ip+ and i/'Z can have zeros 
kj and k/ in the corresponding domains of analyticity: 
det^p+{kj) = 0, kj G C+, and det^Z^{Ki) = 0, k; € C_. 
In virtue of the involution (j3.12p . we obtain k/ — k* and 
equal number N of zeros in both half-planes. The corre- 
sponding RH problem is said to be nonregular, or the RH 
problem with zeros. They are zeros of the RH problem 
that determine soliton solutions of the BEC equations. It 
is seen from Eqs. (|3.6p and p.7p that zeros of ■0+ nullify 
2x2 minors of 4'+- Hence, the rank of ipj^{kj) can be 
equal to one or two. It means in turn that there exist one 
(|lj)) or two (|lj) and \2j)) four-component eigenvectors 
that correspond to zero eigenvalue of ip^ (kj ) : 

ip+{kj)\lj) = for rank?/'+(fcj) = 1, (4.5) 
ip+{kj)\lj) = i/'+(fcj)|2j) = for vankip+ikj) = 2. 

The geometric multiplicity of kj is equal to the dimen- 
sion of the null space of ^+{kj) (1 or 2 in our case). In 
this paper, we only consider the case of zeros kj with its 
geometric multiplicity equal to the algebraic multiplicity 
[which is the order of the zero kj in det ip^{k)]. Note 
that the solution of the RH problem for the general case 
of zeros with unequal geometric and algebraic multiplic- 
ities was elaborated in Ref. [4(j |. 

We will solve the matrix non-regular RH problem 
with zeros fci and fcjf by means of its regularization, 
i.e. by extracting from ipj^. and "02 rational factors 
that are responsible for the appearance of zeros. Hence, 
det^p+{ki) = [and correspondingly det^Z^{ki) = 0]. 
We need a rational matrix function 'E.~^{x, k) which has 
a pole in the point ki. Let us take S~^(a;, k) in the form 



E-\x,k)^l + 



ki - kl 
k — ki 



PM, 



where 



PW= Y. IO(Ai"') 



hn\'m\ 



(4.6) 



/,m— 1 



{m\ = \rfi)'^ due to involution, and r = rank-!/;_|_(fci). P*^'"^ 
is a projector of rank r, (p(''))2 = P^''), and entries of 
the r y. r matrix M are determined by 



V. SOLITON SOLUTIONS 
A. Rank-one soliton 



In the appropriate basis the projector is represented as 
P(i' = diag(l, 0,0,0) or p(2) = diag(l, 1, 0,0). This 
yields 

Therefore, the product V'+(2;, fc)S~^(a;, fc) is regular in k\. 
In the same way, the regularization of V'Z in the point 
k\ is performed by the rational fmiction 



S(.,^) = 1-^P(^), 



(4.7) 



which provides the product E!'0_ to be regular in k\. 
Therefore, the analytic functions are factorized as 

^+(A:) = V^+(fc)S(fc), Vl'W-S-i(fc)^Zi, (4.8) 

with holomorphic functions t\}± which determine the reg- 
ular (without zeros) RH problem: 

V^Z\fc)V^+(fc) = E{k)EG{k)E-'^E-\k), k e Re. 

(4.9) 

For several pairs of zeros {kj,k*), j > I, the regu- 
larization of the RH problem can be performed in the 
same step-by-step manner, with the appropriate defini- 
tion of the eigenvectors within each step. However, for 
practical calculation of A'^-soliton effects it is much more 
convenient to expand a product of rational factors into 
simple fractions, thereby transforming the product-type 
expression into a sum-type one [32l . |40j . 

It is easy to find the coordinate dependence of the 
eigenvectors. Indeed, differentiating (14. 5p in x and in 
t with j = 1 and accounting Eqs. (|2.6p and (|2.7p gives 

(l^i) ^ 10) 



d^\l) = ik,A\l), dt\l)=2ikl\l), 



(4.10) 



with I = 1 for rank one, and 1 — 1,2 for rank two. Hence, 



|/) == eyii>{ikiAx + 2iklAt)\l'-"^), 



(4.11) 



where Ir*^-*) is the coordinate- free four-dimensional vec- 
tor. 

Zeros kj and vectors \lj ) comprise the discrete data 
of the RH problem that determine the soliton content of 
a solution of the BEC equations. The continuous data 
are characterized by the off-block-diagonal parts of the 
matrix G{k) (|4.2[) . k S Re, and are responsible for the 
radiation components. In the following section we con- 
cretize the above relations to obtain one-soliton solutions 
of Eqs. (ESD. 



To obtain the rank-one soliton solution of the BEC 
equations (j2.3p . we consider the single pair fci and fcj of 
zeros and the eigenvector |1). In accordance with Eq. 
(14. lip , the eigenvector takes the form 



^ikix+2iklt ^ikix+2iklt \T 



(5.1) 



where Ua, a = 1, ... ,4 are complex numbers. The RH 
data are purely discrete: N = 1, G{k) = 't, tp± = 1. 
Hence, the solution of the RH problem is given by the 
rational function S in (|4.7p with the projector P'^'^\ The 
reconstruction formula (|4.4p is simplified to 

Q = -2v[K,P^\ (5.2) 

where we set fci = ji + iv and the projector P^^' in 
is explicitly written as 



P 



(1) 



1 



"1 



KP)(KP 



\n,?)V"'"P 



X e~2'^^~'*'(^'-''''*sechz 
Here 

Pab = rianl, z = 2v{x -I- 4^f) -^ p, e 



2p _ \nir + \n2\ 



|n3P + |n4|2- 



Hence, it follows from Eq. (|5.2p that the soliton solution 
is given by 

Q = 2i/H(i)e-2v:r-4»(M'-''')tsechz (5.3) 

with the polarization matrix 



1/2/ niTig nin\ 



n«=i[(K|^-fKn(K|^+|n,|^)] ^^^^^^^^^ 

Note that n2n%^ = nin\ due to the structure of the matrix 
Q in (|2.4p . Besides, the matrix H^^^ obeys automatically 
two conditions: 



dctn(i) = o, |Hii 



(1)|2 



IH 



(1)|2 



2|H«p = l 



Moreover, it is not difhcult to show that the matrix 
H^^) depends only on two essential real parameters. In- 
deed, the rank-one soliton (|5.3p can be represented as 



Q^2v\ ^ „^"" : T^rr. e^'^sechz, (5.4) 
' cos sm e^'^sva ' ^ ' 



where 



cos 9 



l«il 



I'^al 



|ni|2-H|n2|2 |n3|2-H|n4P^ 



X = arg(n3-n4). 



Lp — ~2fj,x — 4(/i ~ ly )t + (pa, 
ipa = arg(ni - rn) = arg(n2 - ns)- 

The soliton amplitude is determined by the parameter 
u, and its velocity is equal to 4/i. The parameters p and 
(j)a give the initial position of the soliton center and its 
initial phase, respectively. The angle 9 determines the 
normalized population of atoms in different spin states, 
while the phase factor e'-^ is responsible for the relative 
phases between the components (p± and (t>o- 

It should be noted for future use that the constant soli- 
ton parameters acquire in general a slow t dependence in 
the presence of perturbation. This results in a modifica- 
tion of the equations for coordinates: 



V 



m = -^(8 / diV(t')^(^')+pW), 



(5.5) 



b{t) = 



2p(,{t)-4: I dt'[ 



^^'it')-ly'{t')] + Mt)■ 



B. Rank- two soliton 



As before, we begin with the pair ki and fc* of zeros, 
but now we have two linearly-independent eigenvectors 

11) = (e~*'^^'^~^''^i*pi, e~*'^^'^~^*'^i*p2, 
|2) = (e-''=i"'-2''=i*gi, e~''=i^-2ifc^^^^ 

ikix+2ik^t„ ikix+2ik^t„ \T 

with Pa and Qa, a — 1,...,4, being complex numbers. 
The rational function S is given by Eq. (|4.7p with the 
rank-two projector P'^' . This projector is written in ac- 
cordance with Eq. (|4.6p as 



2 

E 

l.m—l 



P(2)^ V \m){M-'U{l\^{detMr' 



X (M22|l)(l|-Mi2|l)(2|-M2i|2)(l|+Afn|2)(2|). 
In this case 



M 






z' = 2:^(x -I- 4/xi), and 

Al-blP + b2p, Sl-b3p + b4p, 
A2 = kl|' + 1-721', S2-k3|' + |'?4|', 

A3 = pUi + P2Q2, B3 = pIqs + plq4. 

Introducing the notations (to reproduce literally the re- 
sults of Ref. (3) 

Pl92 -^291 = e''+"^, P3q2 -^293 = /?*, PiQi - PiQi = J*, 
PiQi - PiQi = P4.q2 - P2q4 = a* , (5.7) 



we write explicitly the projector P^"^^ as 



/ Pii P] 



p(2) = 



12 

P12 P22 
P*3 F*4 1 



\P. 



14 -* 24 




where 



= e"^Z-i(ae^-a*I?e~^) 



Pii 
P22 

P12 = -Z-i(a*/?-Ka7*), (5.8) 

-P14 

Pi3 = e"^Z-i(/3e"+7*I?e~"), 

P24 = e"^Z-\je'+P*Ve~'), 

r» = detn^'\ ip = -2px - A{p'^ - iy'^)t + a, 
Z ^detM = l + e^^ + \V\'^e-^\ (5.9) 



n^^) is the polarization matrix II'-^^ 
to the normalization condition [19 



(3 a 
a 7 



1. 



subjected 



(5.10) 



As a result, we immediately find from Eq. (|5.2p with P^^^ 
the rank-two soliton solution of the BEC equations 

M 



Qix,t) ^4:ve"^Z 



n 



(2)g. 



a2Tl'^'^^'' (72T>e' 



, (5.11) 



where CT2 is the Pauli matrix. Notice that the soliton 
solution of the matrix NLS equation was previously ob- 
tained in Ref. [231 by means of the Gelfand-Levitan in- 
tegral equations, while our derivation is purely algebraic. 
The soliton (|5.1ip was also derived by Gerdjikov and co- 
workers via the dressing procedure 41 1 . 

We will distinguish between two featured cases of 
detn(2), namely, detH^^) = and detn^^) ^ 0. These 
cases display different spin properties. Indeed, the spin 
density vector f{x, t) = ti:{Q'^aQ), where ~a is the set of 
the Pauli matrices, is given in general by a spatially odd 
function 



I{x,t) 



(^ 



\V?e- 



a(3 + a(3 -\- a'^ + a^y 
i{a(3 — a(3 + a^ — a^y) 

W-\i? 

with absolute value being of the form 



(5.12) 



l/IH'f 



„2z 



iPPe^'l 



(l-4|X' 



2\l/2 



(5.13) 



Therefore, the total spin vector F = J dxf{x, t) is zero. 
However, for I? = 0, as it follows from Eq. (|5.12p . the 



absolute value of the total spin vector is nonzero, \F\ = 
iiy ^ 0. In accordance with this property, the case V = 
corresponds to the ferromagnetic state, while the case 
2? ^ is usually referred to as a polar state. In fact, a 
true polar state corresponds to the condition |2?| = 1/2, 
when, as it is seen from Eq. (j5.13p . the spin density is 
zero everywhere, not only the total spin [22|. 

It follows from Eq. (|5.1ip that the ferromagnetic state 
has the hyperbolic secant form 



a 0.15 



Qf = 2iyn(2)e''^sech; 



(5.14) 



where entries of the polarization matrix obey the normal- 
ization condition (|5.10p and in addition the constraint 
Pj — a^ = 0. These two condition are sufficient to reduce 
the matrix n*^^^ to the two-parameter form ()5.4p with the 
identifications 




FIG. 1: Profiles of the atom number density function of the 
polar state: \V\ = 1/8 (thick line), i©| = 1/20 (thin line), 



COS 9 = - — p 

(bah 

fa = arg(pi 



bsl 



+ b4P)l/2' 

■Ps) = arg(p4 



X = arg(p3 
-P2)- 



P'l), 




-20 20 



Therefore, the rank-two ferromagnetic soliton is com- 
pletely equivalent to the rank-one soliton (|5.4p . Introduc- 
ing the atom number density n{x,t) and energy density 
e{x,t), 

n{x, t) = tr(QtQ)^ e(x, t) = ctr(Qt Q^ _ Q^QQ^Q), 

as well as their total counterparts Nt — J dxn{x, t) and 
Et = J dxe{x,t), we obtain explicitly the total number 
of atoms and total energy in the ferromagnetic state: 

N^ = Av, E(, = AcN{.{n^~iy^/3). 

In turn, the total number of atoms in the polar state 
and its energy are given by 

NP=8iy, EP.^AcNP{fi^-iyy3). 

The energy difference between both states with equal (ISill and ([7^ with e = 0.1 




-20 20 




FIG. 2: Splitting of the components <l)± and (j)o of a perturbed 
polar soliton and of the atom number density n{x,t). Here 
|/3|2 = 0.7, lap = |7p = 0.1. The perturbation is of the form 



amount of atoms is E^ - E!^ = -{l/16)c{N^)^ < 0. 



Hence, the ferromagnetic state is energetically preferable, 
from the viewpoint of stability, as compared with the po- 
lar state. 

The atom number density of the polar soliton is de- 
scribed by the function 



np{z) 



if 



{e^'+4\V\^ + \V\^e 



2„-2z\ 



(5.15) 



Figure [T] demonstrates typical profiles of the atom num- 
ber density function (|5.15p for different jl^l. The two- 
humped structure becomes more pronounced with de- 
creasing |I?|. Such a state can be treated as a pair of two 
ferromagnetic solitons with antiparallel spins [22l | . Previ- 
ous analysis of stability of multi-humped vector solitons 
for the cubic nonlinearity revealed that they are always 
unstable [42, |43i|. Hence, we can suggest that the most 



likely scenario of the polar soliton evolution under the ac- 
tion of a perturbation would be its splitting into a pair of 
ferromagnetic solitons. Indeed, extensive simulations of 
the perturbed polar soliton behavior demonstrates unam- 
biguously such a splitting. An example of such a behavior 
is depicted in Fig. [21 where we consider a disturbance of 
the intcgrability condition ()2.2p as a perturbation with 
a small parameter e = cq — C2 (see Eq. (j7.3p below for 
a functional form of the perturbation). It is seen that 
all of the components of the polar soliton split under the 
action of the perturbation. 

Let us summarize the main conclusions concerning the 
soliton solutions which will play the key role in study- 
ing soliton perturbations. First, we derived the rank- 
one soliton solution with the hyperbolic secant profile. 



Second, rank-two solutions were obtained and classified 
as ferromagnetic and polar solitons. The polar soliton 
is perturbatively unstable and splits into two rank-two 
ferromagnetic solitons. Third, we proved equivalence of 
the rank-one soliton and rank-two ferromagnetic soliton. 
Therefore, it is sufficient to elaborate a perturbation the- 
ory for the more familiar type of solitons - the rank-one 
soliton (|5.4p . This will be done in the following section. 



VI. PERTURBATION THEORY FOR THE 
BRIGHT SPINOR BEG SOLITON 

In this section we perform a general analysis of the 
perturbed spinor BEC equations 

tdt<l>±+dl^± + 2(|c/.±|2 + 2|0oP)0± 



Then 



2</'^0o = <^^± 



(6.1) 



Here R± and Rq determine a functional form of a per- 
turbation, and e is a small parameter. To distinguish be- 
tween the 'integrable' and 'perturbative' contributions, 
we will assign the symbol 5/5t to the latter. Hence, 



.5Q ^ -^ fO R 



R^ 



R+ Rq 
Rq R— 



In general, a perturbation causes a slow evolution of 
the RH data. Indeed, a perturbation leads to a variation 
SQ of the potential entering the spectral equation (|2.6p . 
and in turn to a variation of the Jost solutions: 

6J±.j; = ik[A, SJ±] + SQJ± + QSJ±. 

Solving this equation gives 

SJ± = J±E ( I Ax'E-^Jg^6QJ±E\ E^K 



As a resuh, we find from Eqs. dSH), ([33]), and (pAO]) a 
variation of the scattering matrix: 



SS ^ '•°° 

- = -^eS^ 



dxE'^i^l^Ri^+ES: 



= -ieT+^ / dxE-^tPZ^R^P^ET^. 

Here S± and T± are the matrices defined in Sect. IIIII 
Notice that they are the analytic solutions ip± that enter 
naturally into this equation. Let us denote 

T±(a,6) = / dxE^^^Pg^RiP^E, (6.2) 

J a 

T±{k) = T±(— oo,oo). 



6S_ 
~5i 



~ieS+T+{k)SZ^ = -ieT-^T_{k)T_ 



The matrices T± are interrelated by means of the matrix 
G entering the RH problem (|4.ip : 



T_(A:)=GT+(/c)G~i. 



(6.3) 



Eventually, variations of the analytic solutions follow 
from Eqs. ([SS]) and ([XTU]) : 



^+ - -ieil>+EH+E-^, 



1„;,-1 



5t 



St 



ieEH^E-'ijjZ 



Here H± are the evolution functionals [30|, |35| that are 
defined in terms of T±, 

H+ = T+(fc)Mi-T+(x,oo), (6.4) 

H_ = MiT_(fc)-T_(x,oo), Ml =diag(l, 1,0,0), 

and contain all essential information about a pertur- 
bation. In particular, the evolution equations for 'ijj± 
gain additional terms caused by the perturbation and 
expressed in terms of H± : 

dtip+ = 2iP[A,i;+] + Vi>+-ie^p+EH+E~\ (6.5) 
dti^Z'^ = 2ik'^[A,i}Z^]-i}Z^V + ieEH^E-^i}Z'^. 

Besides, the evolution equation for the matrix G of the 
RH problem has the form 



dtG = 2iP[A, G] - ieiGH+ - H_G). 



(6.6) 



In fact, this equation gives the evolution of the continuous 
RH data. Note that the involution p.l2p connects i/+ 



with H^ : H^ 



H 



t 



k e Re. 



Now we consider a single rank-one soliton and de- 
rive perturbation-induced evolution equations for the dis- 
crete RH data, i.e. for the zero ki and the eigenvec- 
tor |1). It is more convenient to work with the vector 
\n) = (ni, ri2, ^3, 7^4)"^ which is constant in the absence 
of perturbation and acquires slow t dependence under the 
action of a perturbation. We start from the equation 



ip+{ki)\l) =-0+(fci)exp 



ikix + 2i dtkj A 



b)=0 



which is valid irrespectively of the presence of a pertur- 
bation. Here the integral in the exponent accounts for 
a possible perturbation-induced time dependence of the 
zero fci. Taking the total derivative in t, we obtain 



dt 

+dk 






dtk 



+tl>+(ki)e'^'^'''+^' ^ '^'^^^^dt\p) = 0. 



The first term with dtip+ is given by Eq. (|6.5p which 
contains the evolution functional H^. Recall that the 
evolution functional H^{k) is defined in terms of T+ in 
(|6.4p which in turn depends on V'+^- Hence, the function 
Hj^ is meromorphic in C+ with the simple pole in fci, 
where 1/'+ has zero; 



Just in the same way we obtain evolution equations for 
the parameters tpa and p which are also expressed in 
terms of Ua '■ 



dtfa = 2 



Xn + ^n + (^12^ 



Xl^e'^jtani 



H^{k)^H^-^\k) 



1 



k — ki 



Res[H+{k),ki 



^p+rVa. ^Xii cot 6 + X42e'^) 

gP-»v„ (X*i cot 9 + Xl^e-'^ 



(6.13) 



Here H^^^' stands for the regular part of i?+ in the point 
ki. Following now the method developed in Refs. |30l.[33l. 
[35j . we find that the perturbed evolution of the vector \n) 
is given by 

dt\n) = zee-2^^dtfe?A^|rcg)^^^^^2./dtfc?A|^^^ (g 7) 



dtp^ |{(Xn - X*i) cos2 B + (X22 " ^2*2) sin" & 
+ [{X12 - ^2*1) e''^ - (^1*2 - ^21) e"^^] sin 61 cos 61 



oP+iVc 



{Xj,ie-'^ cos^ 



(6.14) 



Since the left-hand side of Eq. (|6.7p is evidently x- 
independent, we can consider this equation for x — > +00, 
where i/+ has only two non-zero columns [see Eq. (|6.4p ]: 



H+{x 



= (t^}^ ,T^J} A)A) 



(6.8) 

Hence, Eqs. (j6.7p . written in components, take the form 

dtui = ie {XiiTii + Xi2n2) , 

9(712 = ie {X2ini + X22n2) , (6.9) 

dtus = ie {X^iui -I- X32n2) , 

9(714 = ie {Xaitii + ^42772) , 

where for simplicity we use the notation Xab = 
^t£{ki), a,b = 1,2, and X,, = T^l^£ {k,)e-^^ i ''^^l 

for a = 3,4 and b — 1,2. Here T^''^-^ is the regular part 
of T+ in the point fcii 

tV"^'(A:i) = fdxE-^{ki){R{t~P^'^) (6.10) 

+ p(i)i?F(i)+2j/x[A,p(i)i?(l-p(i))]}£;(fci). 

This relation follows from T+(fc) = J dxE-^E-^REE. 
Note that in virtue of the specific symmetry of the matrix 
Q (|2.8p the entries T+32 and T+41 are equal, as well as 
the entries T_|_i4 and T+23. 

Now we can derive the evolution equation for the pa- 
rameters and X entering the polarization matrix of the 
soliton solution (|5.4p . Indeed, these parameters are de- 
fined in terms of 77q which in turn obey Eqs. (|6.9p . Simple 
calculation gives 



2X41 sin 61 cos 61 + Xi2e'^ sin^ 6) 
2X4*1 sin 61 cos 61 + X^^e'"^ sin^ 8) | . 



In fact, Eqs. (|6.1ip - (|6.14p are greatly simplified when 
calculating the functions Xat for a specific perturbation. 
This will be demonstrated in the next Section. 

To derive evolution equation for fci (and hence for the 
soliton amplitude v and velocity /x), we start from the 
equation det^p+{kl) = 0. Taking the total derivative in 
t yields 

dt {detij+{k))^f^^ + {dkdctyj+{k))^^^dtki = 0. 

In accordance with Eqs. (|4.8p and (|4.7p we can write 

/c — fci 



det Tp^{k) 






deti/'+(fc), 



where detip+{ki) 7^ because ip+{k) is a solution of the 
regular RH problem (|4.9p . Accounting now for the rela- 
tion 

9tdet'0+(fc) = — 7etr_ff+ det '0+(fc), 

we eventually obtain a simple evolution equation for the 
zero fci: 



dtki = -etrRes[ff+(fc),fci] 



(6.15) 



■eRes[T+n(fc) + T+22(fc),fci 



dt cos = — 



^P+iVc 



{X^ie-'^ cosO + Xiisml 



eP-"^" {Xl^e''^^ cos 9 + Xl^ sin t 



, (6.11) 



dtX^-^ 



eP+'f" (ATaie-^x - ^426*^ + (tang - cot6l)X4i) 



(X3*ie*x - X*2e-*x + (tan 9 - cot 9)Xl^) 



(6.12) 



Summarizing, Eqs. (HH), ((6lT|) - ((67T4)) . and ([QS]) de- 
termine perturbation-induced evolution of the RH data. 
It should be stressed that these equations are exact be- 
cause we did not yet refer to smallness of e anywhere. 
At the same time, these equations cannot be directly ap- 
plied because T± entering them depend on unknown so- 
lutions V'± of the spectral problem with the perturbed 
potential Q. To proceed further, we develop, owing to 
the smallness of e, the adiabatic approximation of the 
general perturbation theory. 
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VII. ADIABATIC APPROXIMATION 

In the framework of the adiabatic approximation, we 
assume that the perturbed soHton adjusts its shape to 
the unperturbed one at the cost of slow evolution of its 
parameters. Hence, only the discrete RH data are rele- 
vant in this approximation, and we can put ■)/'+ = 1 for 
the solution of the regular RH problem (|4.9p . Therefore, 
-0+ = S. In other words, it is the rational function S that 
completely determines soliton dynamics in the adiabatic 
approximation. In particular, we have 



T, 



dxE-^'~-~^RE.E, E-\k)=E\k). (7.1) 



As an important example, we consider a perturbation 
caused by a small disturbance of the integrability condi- 
tion (12. 2p . In this case we introduce a small parameter as 
e = cq — C2, while the functional form of the perturbations 
_R_j_ has the form 



R±fi = {\h 



»±,0- 



(7.2) 



Inserting the explicit expressions for the soliton compo- 
nents (t>±fl (|5.4p into this equation gives 



R+ = (2i^)3e*('^-x) cos 



R- 

Ro 



|3„i(y+x) , 



^^sech^z, 
'sech'^2, 



(2iy)3e*'^cos6'sin6lsech^z. 



(7.3) 



Matrix elements of T+ which are the main ingredients 
of the evolution equations for the soliton parameters are 
found from Eqs. (|7.ip and (|4.7p . and the projector P^^") 
is calculated by means of the simple formula 



p(i) = 



|1)(1| 
(1|1> 



(l| = |l)t, 



which follows from Eq. (|4.6p . The eigenvector |1) is 
given by Eq. (|5.ip . As a result, matrix elements of the 
projector are as follows {Pta = Pat)' 



P. 



(1) 



1 



5(1) 



1 



e^cos^sechz, P;„' = -e^ ^'^ cos0sin6'sechz, 
pW ^ ie'(¥'-x)cos2 6lsechz, P^V = ie*(^+x) sin^ 6'sech, 



P^W ^ pW ^ _e«Vcos6'sin6'sechz, P^^^ = - 



sin^ 6 sechz, 



.(1) 



.(1) 



P^^' = -e^'^cos^esechz, P^^' = -e""^ sin^ e'sechz, 
P^l^ = ie-^+*>^cos6lsin6'sechz. 



1 lA |2 





FIG. 3: (Color online) Left-hand panel: evolution of the 
perturbed (f)+ component profile obtained numerically. Here 
e — 0.1. Right-hand panel: comparison of the analytically 
predicted frequency shift of the perturbed soliton with that 
obtained numerically from Eqs. (|6.1|l and (|7.2p with e — 0.1. 
In both panels \a\ = |/3| — \'y\ = 0.5. 



their initial values. Finding evolution of the other soli- 
ton parameters demands knowledge of the regular part 
of T_|_(fci). Calculation due to Eq. (|6.10p gives 

-'^ll = -^^12 = -^^21 = -^^22 = 0, 

-'^ai = (2i^^) cxp(-p - iipa + ix) cos^ Q, 
X42 = (2;^2) cxp(-p - iipa - ix) sin^ ^, 
X41 — X-i2 = {2i^^)exp{—p — iipa) cos 9 sin 6. 



Substituting these functions into Eqs. (|6.1ip - (|6.14p . we 
obtain: 

p — const, 9 — const, 

X — const, (pa{t) — ((5q(0) — Aev'^t. 

As a result, within the adiabatic approximation, the only 
manifestation of the perturbation caused by a small devi- 
ation from the integrability condition (j2.2p consists in a 
small shift of the soliton frequency equal to A:tv^ . Hence, 
a ferromagnetic soliton is a pretty robust object against 
a small disturbance of the integrability condition. 

This conclusion has been checked by comparison with 
direct simulations of the perturbed equations (|6.ip . The 
left-hand panel of Fig. [3] demonstrates the evolution of 
the perturbed </)+ component profile. We see a small 
profile distortion. Very little energy radiation is emitted 
to the far field. The same results are valid for the other 
two components. It is seen from the right-hand panel 
that there is a good agreement of the predicted linear 
dependence of the frequency shift on e with that obtained 
numerically. 



VIII. CONCLUSION 



Now we easily obtain from Eqs. (|7.ip and (|7.2p that 

Res[T+n(fc) + T+22(fc),fci]=0. 



Therefore, dtki = in accordance with Eq. (|6.15p . which 
means that the soliton amplitude and velocity preserve 



In the present paper we have developed a perturbation 
theory for bright solitons of the integrable spinor BEC 
model. This model is equivalent to the 2x2 matrix NLS 
equation and is naturally associated with the matrix RH 
problem. We have demonstrated the efhciency of the 
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formalism based on the RH problem, for solving both in- 
tegrable and nearly integrable versions of the spinor BEC 
model. We have obtained the rank-one and rank-two soli- 
ton solutions of the model. Depending on the spin prop- 
erties, the rank-two soliton can be of the ferromagnetic 
type or of the polar type. We have proven that the ferro- 
magnetic soliton is equivalent to the rank-one soliton. As 
regards the polar soliton, its profile is characterized by 
a two-humped structure in a wide region of the soliton 
parameters. We have observed from numerical experi- 
ments that the polar soliton is unstable under the action 
of a perturbation and splits into a pair of ferromagnetic 
solitons. Owing to this fact, the problem to construct a 
perturbation theory for the spinor BEC solitons has been 
reduced to that for the rank-one solitons. 

We have derived perturbation-induced evolution equa- 
tions for the soliton parameters. In the adiabatic ap- 
proximation of the perturbation theory these equations 
have been applied to a practically important case of a 
perturbation caused by a small deviation of the model 
parameters from those in the integrable case. We have 
shown a considerable stability of the ferromagnetic soli- 
ton in the presence of such a perturbation. Namely, the 
soliton preserves its amplitude, velocity and spin proper- 
ties, a small frequency shift being the only manifestation 
of the perturbed environment. At the same time, the 
polar soliton solution of the integrable model has a re- 
strictive area of applicability due to its instability and 
splitting under perturbations. 

Three more points deserve a special comment. First, 



instability of a perturbed polar soliton and its splitting 
into ferromagnetic ones have been observed numerically. 
Analytical study of this phenomenon demands a sepa- 
rate consideration and can be performed, for example, 
by a stability analysis as developed in Ref. J43|. Second, 
we have restricted ourselves to the study of the adia- 
batic approximation of the general perturbation theory. 
Our equations permit us to go beyond this approximation 
and take into account the soliton shape distortion effects. 
However, quantitative characteristics of the first-order ef- 
fects are too small to be verified experimentally, at least 
at present. Examples of practical calculations in the first- 
order approximation can be found in Ref. [33|. Third, 
the formalism developed here for the single perturbed 
soliton can be straightforwardly generalized to the case 
of N weakly interacting solitons arranged into a train- 
like configuration. Analysis of the soliton train dynamics 
by the soliton perturbation theory can be found in Refs. 
[44I . [45J for optical solitons and in Ref. ^] for scalar 
bright BEC solitons. 
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